uji_autokorelasi_spasial

Load Library

library(spdep)
library(tidyverse)
library(readr)
library(readxl)
library(car)
library(sf)
library(stplanr)
library(terra)
library(tmap)
library(leaflet)
library(tseries)
library(dplyr)
library(sfweight) 
library(ggplot2)
library(rgdal)

Load data

dataset <- read_excel("E:/FOLDER KULIAH/STIS/SEMESTER 7/SKRIPSI/dataset_skripsi_revisi.xlsx")

# Check missing data
sapply(dataset, function(x) sum(is.na(dataset)))

##                   No        KodeKecamatan        NamaKecamatan 
##                    0                    0                    0 
##            PEND_LGSG           ntl_median           lst_median 
##                    0                    0                    0 
##           bui_median            co_median         pm2.5_median 
##                    0                    0                    0 
##     elevation_median aridity_index_median         rwi_weighted 
##                    0                    0                    0 
##  jlh_klg_listrik_pln         jlh_fas_pddk         jlh_fas_ksht 
##                    0                    0                    0 
##              jlh_imk  jlh_sentra_industri       jlh_sarana_eko 
##                    0                    0                    0 
##        jlh_tmpt_ibdh   kepadatan_penduduk 
##                    0                    0

str(dataset)

## tibble [572 x 20] (S3: tbl_df/tbl/data.frame)
##  $ No                  : num [1:572] 1 2 3 4 5 6 7 8 9 10 ...
##  $ KodeKecamatan       : num [1:572] 3301010 3301020 3301030 3301040 3301050 ...
##  $ NamaKecamatan       : chr [1:572] "Dayeuhluhur" "Wanareja" "Majenang" "Cimanggu" ...
##  $ PEND_LGSG           : num [1:572] 23.9 16.1 19.8 20.6 25.5 ...
##  $ ntl_median          : num [1:572] 2.06 2.06 1.97 2.55 2.75 ...
##  $ lst_median          : num [1:572] 27.3 28.6 27.1 28.6 28.8 ...
##  $ bui_median          : num [1:572] -1.09 -1.04 -1.04 -1.02 -1.07 ...
##  $ co_median           : num [1:572] 0.0296 0.0307 0.0301 0.032 0.0293 ...
##  $ pm2.5_median        : num [1:572] 40.4 35.1 35.7 33.1 30.3 ...
##  $ elevation_median    : num [1:572] 347 113 214 155 112 52 15 6 3 14 ...
##  $ aridity_index_median: num [1:572] 2.13 1.8 2.04 2.03 1.88 ...
##  $ rwi_weighted        : num [1:572] -0.066 0.4655 0.4854 0.2715 0.0747 ...
##  $ jlh_klg_listrik_pln : num [1:572] 20142 40280 50842 38141 28380 ...
##  $ jlh_fas_pddk        : num [1:572] 119 191 277 180 166 149 136 194 85 208 ...
##  $ jlh_fas_ksht        : num [1:572] 40 75 82 65 44 42 53 48 32 68 ...
##  $ jlh_imk             : num [1:572] 157 357 490 675 162 ...
##  $ jlh_sentra_industri : num [1:572] 0 0 2 2 4 0 0 0 0 0 ...
##  $ jlh_sarana_eko      : num [1:572] 513 839 934 521 503 ...
##  $ jlh_tmpt_ibdh       : num [1:572] 356 585 684 573 492 375 280 421 266 539 ...
##  $ kepadatan_penduduk  : num [1:572] 256 542 842 632 639 ...

# Filter variabel
data_x <- dataset[, -c(1,2,3,4)]
data_y <- dataset[, 4]

# Analisis Deskriptif
summary(data_x)

##    ntl_median       lst_median      bui_median         co_median      
##  Min.   : 1.373   Min.   :22.64   Min.   :-1.20879   Min.   :0.01706  
##  1st Qu.: 2.216   1st Qu.:28.64   1st Qu.:-0.98658   1st Qu.:0.02682  
##  Median : 2.639   Median :30.57   Median :-0.75872   Median :0.02866  
##  Mean   : 3.128   Mean   :30.30   Mean   :-0.77524   Mean   :0.02882  
##  3rd Qu.: 3.182   3rd Qu.:31.93   3rd Qu.:-0.62883   3rd Qu.:0.03111  
##  Max.   :19.035   Max.   :37.88   Max.   : 0.03596   Max.   :0.03821  
##   pm2.5_median   elevation_median aridity_index_median  rwi_weighted    
##  Min.   :17.50   Min.   :   1.0   Min.   :0.8404       Min.   :-0.4047  
##  1st Qu.:32.10   1st Qu.:  25.0   1st Qu.:1.2586       1st Qu.: 0.2410  
##  Median :36.33   Median : 114.5   Median :1.4555       Median : 0.4299  
##  Mean   :35.75   Mean   : 246.5   Mean   :1.5689       Mean   : 0.4561  
##  3rd Qu.:40.20   3rd Qu.: 331.0   3rd Qu.:1.8331       3rd Qu.: 0.6453  
##  Max.   :49.40   Max.   :1816.0   Max.   :2.7325       Max.   : 1.7157  
##  jlh_klg_listrik_pln  jlh_fas_pddk    jlh_fas_ksht       jlh_imk      
##  Min.   : 3516       Min.   : 31.0   Min.   : 11.00   Min.   :  60.0  
##  1st Qu.:15056       1st Qu.:104.0   1st Qu.: 38.00   1st Qu.: 301.8  
##  Median :19732       Median :133.5   Median : 50.00   Median : 558.5  
##  Mean   :21474       Mean   :142.6   Mean   : 56.86   Mean   : 883.8  
##  3rd Qu.:25679       3rd Qu.:171.5   3rd Qu.: 69.00   3rd Qu.:1081.5  
##  Max.   :61913       Max.   :405.0   Max.   :219.00   Max.   :9152.0  
##  jlh_sentra_industri jlh_sarana_eko   jlh_tmpt_ibdh   kepadatan_penduduk
##  Min.   : 0.00       Min.   :  28.0   Min.   : 50.0   Min.   :  117.0   
##  1st Qu.: 0.00       1st Qu.: 556.2   1st Qu.:206.0   1st Qu.:  721.8   
##  Median : 1.00       Median : 789.5   Median :273.5   Median : 1121.5   
##  Mean   : 2.57       Mean   : 854.0   Mean   :294.6   Mean   : 1747.4   
##  3rd Qu.: 3.00       3rd Qu.:1030.0   3rd Qu.:360.5   3rd Qu.: 1736.8   
##  Max.   :54.00       Max.   :3322.0   Max.   :826.0   Max.   :16094.0

summary(data_y)

##    PEND_LGSG     
##  Min.   : 8.309  
##  1st Qu.:20.531  
##  Median :26.146  
##  Mean   :27.604  
##  3rd Qu.:32.614  
##  Max.   :71.308

Uji Moran’s I

# moran.test
W <- read.gal("E:/GEOJSON/QC_MAT.gal", override.id=TRUE)
w <- nb2listw(W, glist=NULL, style="W", zero.policy=NULL)
# write.nb.gal(w$neighbours, "E:/GEOJSON/queen_cont.gal")
Y <- dataset$PEND_LGSG
result_moran <- moran.test(Y, w, randomisation=FALSE, zero.policy=TRUE, alternative="two.sided", 
                           rank = FALSE, na.action=na.fail, spChk=NULL, adjust.n=TRUE)
result_moran

## 
##  Moran I test under normality
## 
## data:  Y  
## weights: w    
## 
## Moran I statistic standard deviate = 7.084, p-value = 1.4e-12
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.1804663118     -0.0017513135      0.0006616349

Simulasi Monte Carlo

set.seed(123)
MC<- moran.mc(Y, w, nsim=599)

# View results (including p-value)
MC

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  Y 
## weights: w  
## number of simulations + 1: 600 
## 
## statistic = 0.18047, observed rank = 600, p-value = 0.001667
## alternative hypothesis: greater

Local Morans’I

oid <- order(Y)
localMI <- localmoran(Y, w)
signifikansi_y <- rep(0,572)
localMI <- cbind(localMI, signifikansi_y)
localMI[,5] <- round(localMI[,5], 4)

for(i in 1:nrow(localMI)){
  if(localMI[i,5] <= 0.05){
    localMI[i,6] <- 1 # SIGNIFIKAN = 1
  } else{
    localMI[i,6] <- 0 # TIDAK SIGNIFIKAN = 0
  }
}


# Jumlah kecamatan yg signifikan
sum(localMI[,6] == 1)

## [1] 73

# Jumlah kecamatan yg tidak signifikan
sum(localMI[,6] != 1)

## [1] 499

View(localMI)

Mapping Local Morans’I P-values

jateng.map <- st_read("E:/GEOJSON/join_new_data.shp",quiet=TRUE)
# MAPPING LOCAL MORAN'S I
listrik <- jateng.map$LISTRIK
listrik_localMI <- cbind(listrik, localMI)

# install.packages("tmap")
library(tmap)
# tm_shape(localMI1) +
#   tm_fill(col = "Ii",
#           style = "pretty",
#           palette = "RdBu",
#           title = "local moran statistics") +
#   tm_borders(alpha = 0.5)

MoranLocalMap = cbind(listrik_localMI, jateng.map)
MoranLocalMap = st_as_sf(MoranLocalMap, crs=st_crs(MoranLocalMap$geometry))
current.mode <- tmap_mode("view") #switch to interactive map

## tmap mode set to interactive viewing

tm_shape(MoranLocalMap) +
  tm_fill(col = "Pr.z....E.Ii..",
          breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
          palette = "-Blues",
          title = "local Moran's I p-values")+
  tm_borders(alpha = 0.5) +
  tm_layout(title="Mapping Local Moran's I P-Value")

Mapping Local Morans’I values

tm_shape(MoranLocalMap) +
  tm_fill(col = "Ii",
          style = "pretty",
          palette = "RdBu",
          title = "local moran statistics (Ii)")+
  tm_borders(alpha = 0.5) +
  tm_layout(title="Mapping Local Moran's I Values")

## Variable(s) "Ii" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.

LISA Cluster MAP 4 quadrant

quadrant <- vector(mode = "numeric", length = nrow(localMI))
DV <- jateng.map$LISTRIK - mean(jateng.map$LISTRIK)
C_mI <- localMI[,1] - mean(localMI[,1])
signif <- 0.05
quadrant[DV >0 & C_mI>0] <- 4      
quadrant[DV <0 & C_mI<0] <- 1      
quadrant[DV <0 & C_mI>0] <- 2
quadrant[DV >0 & C_mI<0] <- 3
quadrant[localMI[,5]>signif] <- 0

MoranLocalMap$quadrant <- quadrant 
colors <- c("#ffffff", "#2c7bb6", "#abd9e9", "#fdae61", "#d7191c")
clusters <- c("insignificant", "low-low", "low-high", "high-low", "high-high")

# knitr::kable(head(MoranLocalMap, n=10))

tm_shape(MoranLocalMap) +
  tm_fill(col = "quadrant", style = "cat", palette = colors[c(sort(unique(quadrant)))+1], 
          labels = clusters[c(sort(unique(quadrant)))+1], popup.vars = c("nmkec")) +
  tm_view() +
  tm_borders(alpha=0.5) +
  tm_layout(title="LISA Map Cluster")

# Jumlah kecamatan yg signifikan
sum(MoranLocalMap$quadrant == 1) 

## [1] 12

sum(MoranLocalMap$quadrant == 2)

## [1] 29

sum(MoranLocalMap$quadrant == 3)

## [1] 13

sum(MoranLocalMap$quadrant == 4)

## [1] 19

# Jumlah kecamatan yg tidak signifikan
sum(MoranLocalMap$quadrant == 0)

## [1] 499

Moran’s Plot

# moran.plot(Y, w)
y_scale <- scale(Y)
y_scale <- as.vector(y_scale)
mp <- moran.plot(Y, w, labels=as.character(dataset$NamaKecamatan), pch=16)

# moran.plot(as.vector(scale(Y)), w,
#            labels=as.character(dataset$NamaKecamatan), xlim=c(5, 75), ylim=c(10,60), pch=16)
if (require(ggplot2, quietly=TRUE)) {
  xname <- attr(mp, "xname")
  ggplot(mp, aes(x=x, y=wx)) + geom_point(shape=1) + 
    geom_smooth(formula=y ~ x, method="lm") + 
    geom_hline(yintercept=mean(mp$wx), lty=2) + 
    geom_vline(xintercept=mean(mp$x), lty=2) + theme_minimal() + 
    geom_point(data=mp[mp$is_inf,], aes(x=x, y=wx), shape=9) +
    geom_text(data=mp[mp$is_inf,], aes(x=x, y=wx, label=labels, vjust=1.5)) +
    xlab(xname) + ylab(paste0("Spatially lagged ", xname))
}

Membuat cluster map LISA

# shape <- st_read(system.file("shape/nc.shp", package="sf")) # included with sf package
queen <- poly2nb(as(jateng.map, "Spatial"), queen = T)
# summary(queen)

# calucualte the lisa groups
shape_lisa <- jateng.map %>% 
  mutate(nb = st_neighbors(geometry),
         wts = st_weights(queen),
         lag_listrik = st_lag(LISTRIK, queen, wts),
         lisa = categorize_lisa(LISTRIK, lag_listrik))

## Warning: `st_neighbors()` was deprecated in sfweight 0.0.0.9002.
## i Please use `st_contiguity()` instead.

# report results
ggplot(data = shape_lisa) +
  geom_sf(aes(fill = lisa))

Konektivitas

library(rgdal)
jatengmap <- readOGR("E:/GEOJSON/join_new_data.shp")

## OGR data source with driver: ESRI Shapefile 
## Source: "E:\GEOJSON\join_new_data.shp", layer: "join_new_data"
## with 572 features
## It has 32 fields
## Integer64 fields read as strings:  No KodeKec

plot(jatengmap, border = 'black')
plot(queen, coordinates(jatengmap), pch = 19, cex = 0.6, add = TRUE, col = "red")

rswm_q <- nb2listw(queen, zero.policy = TRUE)

summary(rswm_q)

## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 572 
## Number of nonzero links: 3136 
## Percentage nonzero weights: 0.9584821 
## Average number of links: 5.482517 
## Link number distribution:
## 
##   2   3   4   5   6   7   8   9  10  11 
##   7  41 102 142 146  90  30   6   7   1 
## 7 least connected regions:
## 83 260 295 439 446 454 548 with 2 links
## 1 most connected region:
## 288 with 11 links
## 
## Weights style: W 
## Weights constants summary:
##     n     nn  S0       S1       S2
## W 572 327184 572 218.5639 2336.428